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In lattice QCD, the Maximum Entropy Method can be used to reconstruct spectral functions from 
euclidean correlators obtained in numerical simulations. We show that at finite temperature the 
most commonly used algorithm, employing Bryan's method, is inherently unstable at small energies 
and give a modification that avoids this. We demonstrate this approach using the vector current- 
current correlator obtained in quenched QCD at finite temperature. Our first results indicate a 
small electrical conductivity above the deconfmement transition. 

PACS numbers: 12.38.Gc Lattice QCD calculations, 12.38.Mh Quark-gluon plasma 



In the deconfined, high-temperature phase of Quantum 
Chromodynamics, the behaviour of spectral functions of 
conserved currents at small energies is of intrinsic inter- 
est due to its relation with transport properties of the 
quark-gluon plasma (QGP). According to the Kubo for- 
mulas [l|, transport coefficients, such as the shear and 
bulk viscosities and the electrical conductivity, are pro- 
portional to the slope of appropriate spectral functions at 
vanishing energy. The success of e.g. ideal hydrodynam- 
ics in heavy ion phenomenology, assuming vanishing vis- 
cosities and requiring early thermalization 2] , has lead to 
the notion 0] that the QGP created in relativistic heavy 
ion collisions at RHIC is strongly coupled and that the 
ratio of shear viscosity to entropy density in this sQGP 
may be close to the conjectured lower bound Q reached 
in thermal field theories that admit a gravity dual [B| . 

In order to put these ideas on firm footing, it is im- 
portant to have a first-principle calculation of transport 
coefficients in the strongly coupled regime of hot QCD. 
As is well known 0, 0], a nonperturbative calculation 
using lattice QCD is difficult due to the necessity to 
perform an analytic continuation from imaginary to real 
time. The most common approach used to obtain spec- 
tral functions from euclidean correlators is the Maximum 
Entropy Method Q, employing Bryan's algorithm 
Our first aim in this Letter is to discuss this method in 
some detail and point out a source of numerical instabil- 
ities present in most finite-temperature studies available 
to date. We show how the algorithm can be modified to 
avoid this problem. Our second aim is to apply the new 
method to the vector current-current correlator, obtained 
in quenched lattice QCD simulations at finite tempera- 
ture, using staggered quarks. We study the behaviour at 
small energies and argue that it allows us to extract a 
value for electrical conductivity in the strongly coupled 
regime above the deconfmement transition. 

Maximum Entropy Method - The relation between the 
euclidean correlator G(r) = J d 3 x (J(r, x) Jt (0,0)) at 



zero momentum and the corresponding spectral function 
p(ui) reads 
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G(r)= / -ifKr)p( W ), 
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where the kernel is given by 



K(u,t) = 



cosh[w(r- 1/2T)] 
sinh(w/2T) 



(1) 



(2) 



We consider (local) meson operators of the form 
J(t, x) = q(r, x)rg(r, x), where T depends on the chan- 
nel under consideration. The temperature T is related to 
the euclidean temporal extent N T by 1 /T — aN T , where a 
is the (temporal) lattice spacing. The difficulty in invert- 
ing relation fl} is due to the fact that G(r) is obtained 
numerically at a discrete set of points r, = T m - m + (i — l)a 
(i = 1, . . . , N), where the number of data points N is 
typically O(10), whereas p(u>) is in principle a continuous 
function of uj. Simple properties of the kernel and spec- 
tral functions allow us to cutoff the to integral at w max . 
The resulting finite interval is discretized as w„ = nAw 
(n = 1, . . . , N u ), where N w is typically O(10 3 ), making a 
simple inversion ill-defined (below we suppress the index 
n in w„ where possible). 

Using the ideas of Bayesian probability theory, one 
may construct the most probable spectral function by 
maximizing the conditional probability P[p\DH], where 
D indicates the data and H some additional prior knowl- 
edge. In the Maximum Entropy Method (MEM), the 
prior knowledge is encoded in an entropy term, 



s 



du> 
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p{ijj) — m{ijj) — p(u>) In 



p(u>) 



m{uj) 



(3) 



where m(u>) is the so-called default model, containing the 
additional information. An often used default model is 
m(u>) — mow 2 , with too a (channel-dependent) constant. 
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This choice is motivated by the large-u behaviour of me- 
son spectral functions in the continuum theory, which is 
accessible in perturbation theory. The conditional prob- 
ability to be extremized reads 



xS 



(4) 



where % 2 is the standard likelihood function and a is a 
parameter balancing the relative importance of the data 
and the prior knowledge. Since p(u>) is nonnegative for 
positive id, it is written as 



p(w) = m(w)exp/(w), 



(5) 



where f(uj) is to be determined. 

Bryan's approach - To make the problem well-defined, 
the arbitrary function f(u>) has to be reduced to one con- 
taining at most N parameters, i.e. it has to be restricted 
to an N dimensional subspace. In Bryan's algorithm the 
subspace is determined using a singular value decompo- 
sition (SVD). After discretizing the u> integral, the ker- 
nel K(w n ,Ti) is viewed as a N u x N matrix. The SVD 
theorem states that it can be written as K = UWV T , 
with U an N w x N matrix satisfying U T U — ^nxn, 
W = diag(u>i, . . . , wn) with wi > . . . > wn > 0, and V 
an N x N matrix satisfying VV T = V T V = l N xN- The 
dimension N s of the subspace is determined by the singu- 
lar values Wi. Due to the reflection symmetry K{u>, l/T— 
t) = K(uj,t), one finds that N s < N < N T /2 + 1. Dis- 
carding the contact term at r = 0, we use below the max- 
imal range N s = N = N T /2, such that all data points at 
1 < T; L j a < N T /2 are included. In Bryan's algorithm the 
subspace is defined as the space spanned by the column 
vectors of U, i.e. the N vectors Uj (i = 1, . . . , N) with 
elements Uj(w n ) = U n i. Since U is orthogonal, these vec- 
tors are linearly independent and normalized, with the 
inner product (ui\uj) = J2 n =i u i{ UJ n)Uj{^n) = It 
follows from the extremum condition that /(cj) can be 
parametrized as 



N 

/m = y^cjUi(uj). 

i=l 



(6) 



This reduces the problem to the determination of N co- 
efficients Ci, as desired. 

The Ui functions should therefore be regarded as basis 
functions and it is interesting to study how they behave. 
In Fig.[T]we show the first four basis functions determined 
using the SVD of K(oj n ,Ti) for the case that N u = 1000 
and N = N T /2 = 12. The cutoff is au max = 5. Details of 
the basis functions depend on the choice of N u and N, 
but we find that the i'th function crosses zero i times. 
The inset shows a blowup of the small energy behaviour. 
It can be seen that all basis functions seem to diverge in 
the small energy limit (although they are still normaliz- 
able). This is because the kernel itself diverges for small 
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FIG. 1: First four basis functions Ui(ui) as a function of aid for 
aw max = 5, Nui = 1000, N T = 24, using the standard kernel. 
The inset shows a blow-up of the small energy region. 



id, since 



lim Kkd, t) 



2T 



Ld 
T 



- tT(1 - tT) 
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(7) 

Due to this divergence, it is not possible to include the 
point at Ld = (note that the limits T — > and id — * 
do not commute). 

In the MEM analysis of lattice correlators at finite tem- 
perature, it is often found (see e.g. Ref. [13] and references 
therein) that the small energy behaviour is unstable: the 
value of the spectral function at the smallest nonzero u> 
value is inconsistent with values at larger energies and 
the u — > behaviour depends on Alo, indicating that 
it is an artefact of the method. Moreover MEM does 
not always converge. Unstable features in the algorithm 
at small energies prevent of course insight in transport 
properties of the QGP. 

Modification of Bryan's algorithm - As indicated 
above, the divergence of the kernel at small energies can 
lead to a numerically unstable algorithm. Fortunately, 
this can easily be avoided by writing 



K(Ld,T) = —K(ld,T) 



2T 

pM = — pM 

UJ 



(8) 



such that K(Ld,r)p(Ld) — K (a;,r)/o(a;). We may now re- 
peat the SVD for K(uj n ,Ti). We find that both kernels 
give the same identification of the dimension of the sub- 
space. The first four new basis functions Tii are shown 
in Fig. O The inset shows again a blow-up. Since 
K(0,t) = 1, the basis functions take a finite value at 
small u and Tii(0) is well-defined for all i. From now 
on we include the point at u = in the analysis. The 
redefined spectral function is parametrized as [181 ] 



N 



p{Ld) = m(u) expy^^gui(Ld), 



(9) 



3 



0.3 



redefined kernel 
— u/co) 




60 



FIG. 2: As in Fig. [T] using the redefined kernel K(lu,t). 
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FIG. 3: Vector spectral functions p{u)/uoT as a function of 
lo/T. We used N u = 1000, aw max = 5, 6 = 1. 



with the default model rn(uj) ~ m(uj)/uj. In order to 
make contact with previous results, rfi(uj) ~ u> at large w. 
For small o>, we note that we reconstruct p ~ p/^j such 
that the intercept at uj — is proportional to the appro- 
priate transport coefficient. Specifically, for the electrical 
conductivity we find that 



— = lim 



p(w) 75(0) 



o 6loT 12T 2 '- 



(10) 



where in this case p is in the vector channel (r = 7fc, 
summed over k = 1,2, 3). To allow for a nonzero inter- 
cept, we find from Eq. ^ that m(0) should be finite and 
nonzero. We use therefore the following default model 



a 2 m(uj) — mo(b + aui), 



(11) 



where b < 1 is a parameter that can be varied to assess 
the default-model dependence of spectral functions in the 
low-energy regime. 

Lattice QCD - We now apply the modified algorithm to 
the euclidean correlator in the vector channel, obtained 
in quenched lattice QCD simulations at finite tempera- 
ture. Lattice details are given in Table|TJ The two lattices 
above T c differ only in temperature, while the lattice be- 
low T c is coarser but has N T = 24, in common with the 
hot one. We use light staggered quarks with a bare mass 
of am = 0.01. As is well known, the pure gauge action 
has a global Z3 symmetry, which is broken by the pres- 
ence of fermions in full QCD. To incorporate this, we 





P 


a' 1 (GeV) 




T/T c 


# conf 


cold 


6.5 


4.04 


48 3 x 24 


0.62 


100 


hot 


7.192 


9.72 


64 3 x 24 


1.5 


100 


very hot 


7.192 


9.72 


64 3 x 16 


2.25 


50 



TABLE I: Lattice parameters. Estimates for the lattice spac- 
ing and temperature are taken from Ref. [Toj ] . 



multiply the link variables by an element of Z3 in the 
calculation of the quark propagators so that the phase 
of the Polyakov loop is approximately real. Chiral sym- 
metry restoration at T > T c is then clearly visible: the 
pseudoscalar and scalar correlators coincide [12]. 

For staggered fermions the MEM analysis is compli- 
cated due to the mixing of two signals in the correlator. 
The equivalent of relation ([1]) reads 



G stag (r) 



^K(t,cj) pU.)-(-l) 



7»J5 



p{w) 



(12) 

where p is related to p by changing r to T — 7475T. In 
practice we perform an independent MEM analysis on 
the even and odd time slices, obtaining p ovcn — p — p and 
Podd = P + P and combine these to get p. Let us remark 
here that while the original formulation, using K(ui,t), 
failed to converge in some cases, we found that the new 
method worked successfully (we have studied the vector 
and pseudoscalar channels [13j). In cases where both 
methods worked we found comparable results for spectral 
functions in most of the energy range but deviations in 
the low-energy region, as expected. 

In Fig. [3] we show the vector spectral functions p(ui) 
normalized by wT as a function ofu/T for the three tem- 
peratures. Our treatment of a follows Ref. Q. The large 
uj behaviour at auj > 3 is determined by the continuum 
default model, p{uj)/uj ~ m^u (recall that ui/T — N T auj). 
In most MEM studies we carried out we find "bumps" at 
1 < auj < 2: within uncertainties these do not depend on 
the channel under consideration, quark mass or external 
We interpret them as lattice arte- 



12, LL 



momentum 

facts present due to the difference between continuum 
and lattice fermion dispersion relation (see Ref. llj for 
an analysis of free staggered quarks). The peak structure 
at Lo/T ~ 9 on the cold lattice is the rho- meson, which is 
more clearly visible in plots of p{uj)/u} 2 vs. ui/T [l^ . The 
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FIG. 4: Default model dependence of p(u)/ujT for JV T = 
24 (hot) and 16 (very hot) in the low-energy region. We 
show results for iV^ = 1000, 2000 and b = 1.0, 0.5, 0.1 at fixed 

fik-'max — 5. 



rho-peak is less pronounced on the hot lattice and seems 
to have vanished on the very hot lattice. Error bars are 
explained in Ref. Q. 

Above T c a nonzero intercept at to — can be seen. 
In the cold phase the intercept is an order of magnitude 
smaller. A blow-up of p(uj)/u>T is shown in Fig.SJ In the 
high-temperature weakly-coupled theory, the transport 
contribution is located at ui/T ~ <? 4 (T) <C 1, where the 
spectral function behaves as p(ui)/u>T ~ l/g 4 (T) 3> 1. 
As a result the euclidean correlator is insensitive to the 
details of the spectral function Here, by contrast, 
we find a clear signal for nonzero spectral weight varying 
smoothly in the range < ui/T < 4. To assess this, we 
have varied N u , aw max and the parameter b in Eq. (|11[) . 
In the case that we take 6 = (therefore disallowing a 
nonzero intercept), we find that the MEM routine does 
not converge. Results with N w = 2000 cannot be dis- 
tinguished from those obtained with N u — 1000. These 
findings suggest that nonzero spectral weight and a finite 
intercept are robust outcomes. Dividing the intercept at 
oj = by 6 yields the electrical conductivity. We find 
er/T = 0.4 ± 0.1, where the error is an indication of the 



uncertainty in the MEM reconstruction [14j [19(. Statis- 
tical uncertainty and the effect of operator normalization 
[lit are expected to be smaller. 

Outlook - We identified a source for numerical instabil- 
ities in the standard formulation of MEM with Bryan's 
approach at finite temperature and showed how this can 
be avoided. We then applied this new method to the vec- 
tor current-current correlator obtained in quenched lat- 
tice QCD at finite temperature. We found a robust signal 
for a small but nonzero value of the electrical conductiv- 
ity above the deconfinement transition, in line with the 
notion of the sQGP. Directions for the future are many. It 



is important to repeat the analysis using Wilson fermions 
and the conserved point split current, and to include the 
shear and bulk viscosity as well. The effect of dynam- 
ical quarks can be investi gat ed in two-flavour QCD on 
highly anisotropic lattices [l(| • We emphasize once more 
that the new formulation is essential for analyzing the 
low-energy region. 
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